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We investigate the ability of frame-invariant amplitude equations [G. H. Gunaratne, Q. Ouyang, 
and H. Swinney, Phys. Rev. E 50, 2802 (1994)] to describe quantitatively the evolution of polycrys- 
talline microstructures and we extend this approach to include the interaction between composition 
and stress. Validations for elemental materials include studies of the Asaro-Tiller-Grinfeld morpho- 
logical instability of a stressed crystal surface, polycrystalline growth from the melt, grain boundary 
energies over a wide range of misorientation, and grain boundary motion coupled to shear defor- 
mation. Amplitude equations with accelerated strain relaxation in the solid are shown to model 
accurately the Asaro-Tiller-Grinfeld instability. Polycrystalline growth is also well described. How- 
ever, the survey of grain boundary energies shows that the approach is only valid for a restricted 
range of misorientations as a direct consequence of an amplitude expansion. This range covers 
approximately half the complete range allowed by crystal symmetry for some fixed reference set of 
density waves used in the expansion. Over this range, coupled motion to shear is well described by 
known geometrical rules and a transition from coupling to sliding motion is also reproduced. Am- 
plitude equations for alloys are derived phenomenologically in a Ginzburg-Landau spirit. Vegard's 
law is shown to be naturally described by seeking a gauge invariant form of those equations under 
a transformation that corresponds to a lattice expansion and deviations from Vegard's law can be 
easily incorporated. Those equations realistically describe the dilute alloy limit and have the same 
flexibility as conventional phase-field models for incorporating arbitrary free-energy/composition 
curves. As a test of this approach, we recover known analytical expressions for open-system elastic 
constants [F. C. Larche and J. W. Cahn, Acta metall. 33, 331 (1985)]. 
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I. INTRODUCTION AND SUMMARY 

Rapid advances in phase-field modeling over the 
last two decades have greatly enhanced our ability to 
model a wide range of complex interfacial patterns in 
materials^^^. In mature applications such as dendritic 
solidification, it has been possible to bridge successfully 
atomistic and continuum scales by linking molecular dy- 
namics and phase field simulations'^^S. This bridge has 
relied on the combination of thin-interface asymptotic 
analyses of phase-field modelsiSrJ^ to simulate interface 
dynamics on experimentally relevant length and time 
scales, and of new atomistic simulation methods to pre- 
dict some key parameters for those problems such as the 
anisotropy of the crystal-melt interfaceS.. 

Despite this progress, simulating the evolution of poly- 
crystalline patterns has remained challenging. Those 
patterns have been modeled using multiple phase 
fields^!^, each representing a different crystal orientation, 
or by introducing a scalar order parameter that repre- 
sents the local crystal orientatio n^^d"^ . At a more mi- 
croscopic level, the phase-field-crystal (PFC) approach 
has emerged as an attractive alternative^^"^^. The PFC 
mode l^^d^ is a reformulation of the Swift-Hohenberg 
(SH) model of pattern formation^^ with conserved dy- 
namics and can also be motivated as a simplified ver- 
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sion of classical density functional theory (DFT)i^. By 
modeling directly the crystal density field, it provides 
a simple frame invariant description of polycrystals and 
it naturally incorporates defects and elastic interactions. 
However, resolving the density field on the scale of the 
lattice spacing limits the system sizes that can be stud- 
ied. Therefore, finding ways to "coarse grain" spatially 
the PFC model while retaining its advantages is highly 
desirable, both for computational purposes and to gain 
analytical insights into the model properties. The am- 
plitude equation approach^L^, widely used in a pattern 
formation context;^"— has been revived recently as a 
method to coarse-grain the PFC model^"— . Amplitude 
equations have also been used to study the properties of 
crystal-melt interfaces^. 



A. Frame-invariant amplitude equations 

The amplitude equation approach was pioneered by 
Newell- Whitehead-Segel (NWS)2ii2^ to model stripe pat- 
terns in Rayleigh-Benard convection. While NWS de- 
rived an amplitude equation from a multiple scale anal- 
ysis of the Navier-Stokes equations, the same equation 
can be derived from the SH model where the pattern is 
described by some scalar field For stripes in one di- 
mension, the amplitude equation approach exploits the 
fact that the pattern is slowly modulated in space close 
to the onset of instability, with the distance from onset 
measured by a dimensionless parameter e. Hence ^' can 
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be written as a sum of plane waves 

* = u{X)e"^'>'^ +u*{X)e- 



(1) 



where the complex amplitude u{X) varies spatially on the 
slow scale "X — e^/^qox" as opposed to the original fast 
scale qox, which are both defined here to be dimension- 
less. In addition, qq ~ 271 /a where a, the wavelength of 
the stripe pattern, is the analog of the "lattice spacing" . 
In this framework, coarse graining consists of deriving 
an amplitude equation, i.e. an evolution equation for m, 
which can be solved on the slow scale X. The NWS am- 
plitude equation, which is derived formally from a multi- 
ple scale analysis that exploits the smallness of e, obeys 
the gradient dynamics ii ^ —6F/Su* with the Lyapunov 
functional ^ J dx [\dxu\'^ + fb{u)] and a bulk energy 
term fb{u) that wiU be specified later. 

As a consequence of this coarse graining, the NWS 
amplitude equation is not frame invariant (in contrast 
to the dynamical equation for 'I') because of the fixed 
choice of reference axis for the underlying plane waves. 
To overcome this limitation, Gunaratne, Ouyang, and 
Swinney (GOS)'^^ have derived a more general frame- 
invariant form of the NWS equation in the application 
of this approach to two-dimensional hexagonal patterns. 
For a "crystalline" pattern, 4' can be written as 



N 



vl/(r) ^ u^^^(f) exp{ik^'^ ■ r), 



(2) 



where the sum of plane waves is taken here over all princi- 
pal reciprocal lattice vectors fc'-'-' (with = 6 for hexago- 
nal ordering) and the complex amplitudes have the prop- 
erty that u^^^ = (w^'))* if k'^^^ — —k^^\ which ensures 
that * is real, and Ifc^-''] = go for ah j. The COS ampli- 
tude equations (written here in terms of the fast spatial 
variable f — xx + yij -\- zz, using the hat for normalized 
vectors) introduce the "box operator" 



Uj = ■ V - — v^ 

2^0 



(3) 



which makes those equations rotationally covariant. 
They have the variational form (w*^^^ and its complex con- 
jugate u^J^* are considered as independent fields) 



aw^j) = -r 



where the functional 



T= dr 
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and the bulk free-energy density fb is generally a sum of 
products of amplitudes for which the sum of reciprocal 
lattice vectors form closed polygons up to quartic terms. 
More rigorous derivations of frame-invariant amplitude 



equations using a renormalization group (RG) framework 
have been given for the SH equation^i^'*'* and, more re- 
cently, the PFC model^"— . All those analyses recover 
essentially the same form of the GOS amplitude equa- 
tions with the box operator. The RG analysis of the 
PFC model yields amplitude equations with higher order 
spatial derivatives^^i^ than GOS. However, as clarified 
recently by Chan and Goldenfeld-? , the RG amplitude 
equations for PFC can be derived from the same free- 
energy functional ^ with a form that conserves the vol- 
ume integral of 5*. In the present notation, this conserved 
dynamics is given by 



dtU<-J^ = -r (1 - 2iqQ'D,] 



5T 



(6) 



When expressed in terms of the slow variable R = 
e^^^qor, the operator 1 — 2iqQ^Dj « 1 in the small e 
limit of those equations. Hence the non-conserved and 
conserved dynamics defined by Eqs. (|3]) and (|n|) are es- 
sentially identical in this limit. This limit is physically 
relevant since e has been shown to be small in fits of the 
amplitude equations to actual materials (e.g. pure Fe, 
see Ref. l42h . Therefore, in the present work, we use the 
non-conserved form (|4]) that is both accurate and com- 
putationally more efficient in this small e limit. 

Simulations to date support the feasibility of using an 
amplitude equation approach to simulate polycrystalline 
microstructural evolutional^. Furthermore, they have 
demonstrated the possibility of significant computational 
cost saving by using adaptive meshing algorithms^. De- 
spite this progress, the quantitative validity of this ap- 
proach for modeling polycrystalline patterns is still not 
fully explored. In addition, for materials application, it 
would be desirable to extend this approach to alloys. The 
dual goal of the present work is to explore in more depth 
both analytically and numerically the quantitative valid- 
ity of frame-invariant forms of the amplitude equations 
for simulating polycrystalline pattern evolution and to 
extend this approach to binary alloys. 



B. Extension to alloys 

One possible approach to model binary alloys is to de- 
rive amplitude equations directly from a PFC model with 
two coupled conserved fields^. This approach, which was 
pursued recently by Elder et ali^, yields qualitatively 
similar eutectic phase-diagrams as earlier conventional 
phase-field models of two-phase growth--'--, and has 
been shown to model complex patterns with defects and 
elasticity^. One limitation is that it does not realistically 
describe the dilute alloy limit and more generally lacks 
the fiexibility of the conventional phase-field approach to 
model arbitrary liquid and solid free-energy/composition 
curves^. 

An alternate approach, which is pursued here, is 
to write down amplitude equations phenomenologi- 
cally based on physical considerations in the spirit of 
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Ginzburg-Landau theory. For pure materials, this ap- 
proach was developed by Shih et al.'^^ in the framework 
of DFT to model body-centered-cubic (bee) /liquid inter- 
faces. Even though there has been subsequent attempts 
to derive phase field models from DFT^Ai^, the work 
of Shih et alA^ was among the first to derive an analyti- 
cal form for the "double- well potential" of the phase field 
model and the surface energy in the isotropic limit where 
all density waves have equal amplitudes. In this limit, 
crystalline order is described by a single scalar variable 
directly analogous to the phase field. 

The amplitude equations of Shih et alA^ were recently 
revised by Wu et a/.— with improved predictions bench- 
marked against molecular dynamics simulations for pure 
Fe. In this revision, the coefficient of the gradient square 
term in ([5]) (and quadratic terms in /;,) were related di- 
rectly to liquid structure factor properties by a small gra- 
dient expansion similar to the one used by Haymet and 
Oxtoby in a more complete DFT study of crystal-melt 
interfaces^2ii^. This expansion yields a free-energy func- 
tional of the form ([5]) with Oj « k^^'^ ■ V, where this 
truncation is accurate for the purpose of computing the 
solid-liquid interfacial free-energy and its anisotropy, as 
shown here for the same Fe parameters as in Ref . \^ The 
coefficients of higher order nonlinearities (cubic or quar- 
tic) in fh were obtained using the same ansatz as Shih et 
alM: This ansatz holds that all products of amplitudes 
corresponding to polygons with the same number of sides 
(three or four) have equal weight Wu and Karma''^ 
have shown that this ansatz yields different quartic non- 
linearities than in the amplitude equations derived from 
the PFC model, but that those differences do not alter 
significantly solid-liquid interface properties. 

In extending this approach to alloys, two distinct ef- 
fects of solute addition need to be considered. The first is 
the coupling between composition and crystalline order. 
This coupling can be introduced phenomenologically by 
defining a real scalar function tf>{{u^J^}) of the complex 
amplitudes, which varies between zero in the liquid to one 
in the completely ordered solid. Several possible choices 
of functions that accomplish this goal will be specified 
later in the paper. This crystalline order parameter can 
then be used to interpolate between the thermodynamic 
properties of the solid and liquid phases, as in the conven- 
tional phase- field approach^^, by writing the free-energy 
density due to solute addition in the form 



(7) 

where fs{c,T) and fi{c,T)) are the solid and liquid free- 
energy /composition curves, respectively, T is the tem- 
perature, and c is defined here to be the mole fraction 
of B in A. To construct an alloy model, fc needs to be 
added inside the square brackets in Eq. ^ , while making 
at the same time the bulk free-energy density of the pure 
material, ft, dependent on temperature. 

The second is the coupling between composition and 
stress. A detailed thermodynamic treatment of the cou- 



pling between composition and stress has been given by 
Larche and Cahn (see Ref. [521 with earlier references 
therein). Solute addition generally modifies the equilib- 
rium lattice constant of an alloy. In its simplest form, the 
relationship between the lattice constant and concentra- 
tion at fixed temperature is described by Vegard's law^, 
which is a linear relationship of the form 



ao(l 



aic, 



(8) 



where qq and ai are the lattice constants of pure A and 
pure B, respectively. While this empirical law holds ap- 
proximately for ionic crystals, metallic alloys show sig- 
nificant deviations from this law. The origin of those de- 
viations has been widely studied theoretically including 
in the context of classical DFT—. 

To see how to incorporate the coupling between com- 
position and stress in the amplitude equations, consider 
a simple one-dimensional crystal (stripe pattern) repre- 



sented by 



ue 



iqox 



u*e ^10^ . Since considerations 



of frame invariance are irrelevant in this case, the pure 
limit of this problem is described by the NWS amplitude 
equation with a term in the free-energy density. 

A change of lattice constant (wavelength of 4*) can be 
generally represented by the transformation u — >■ ue^P^^\ 
Assuming that p{x) is slowly varying spatially on the lat- 
tice scale, this transformation changes locally the lattice 
constant to a = 27r/(go + dp/dx), which can be approx- 
imated by a w ao(l — QQ^dp/dx), where ao = 27r/go- 
This approximation is valid as long as the change of lat- 
tice constant is small qQ^\dp/dx\ ^ 1, where dp/dx < 
corresponds to a lattice expansion. 

The transformation u — > ue*'''^-' is directly analogous 
to the gauge transformation of the quantum mechanical 
wavefunction generally considered in deriving a gauge in- 
variant form of the Schroedinger equation for a charged 
particle in an electromagnetic field. Clearly, like the stan- 
dard Schroedinger equation, the amplitude equation with 
a term ~ d^u, derived from the "kinetic" part IcJ^wp 
of the free-energy density, is not invariant under this 
transformation. However, this analogy suggests that a 
gauge-invariant form can be obtained by the substitu- 
tion dx —?' dx + ioic, which transforms the kinetic part 
into \{dx+icec)u\'^ , where a is a coupling constant. The 
new amplitude equation for u (considered independently 
from the governing equation for c) is now invariant under 
the transformations u — ?> ue*''^^-' and c — > c — a~^dp/dx, 
which yields essentially Vegard's law for the choice of 
coupling constant a = qoio-i — ao)/<^o- 

The generalization to more realistic two- and three- 
dimensional crystal structures is immediate since a 
change of lattice parameter corresponds to a change of 
magnitude of each reciprocal lattice vector fc(j), and is 



obtained by the substitution Dj Dj 



Higher 



order spatial derivatives in the box operator give a neg- 
ligible contribution since p{x) is slowly varying on the 
lattice scale. In addition, deviations from Vegard's law 
can be modeled by the more general transformation 
Dj — ^ Dj + iqQ^h{c). This transformation yields the 
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equilibrium lattice constant a ^ ao{l + h{c)). Therefore 
it can describe an arbitrary relationship between the lat- 
tice constant and composition since h(c) can be chosen 
to be an arbitrary function of c. 

As an analytical validation of our approach, we derive 
expressions for the modification of the so-called "open- 
systems" elastic moduli, which agree with the expressions 
derived by Larche and CahrtS. The existence of those 
moduli stems physically from the fact that the crystal 
strain is generally a function of both composition and 
stress. In contrast, the composition for an "open system" 
at fixed chemical potential fj,, in contact with a reservoir 
of solute, is only a function of stress. Hence, the strain 
at fixed n can be expressed as a function of stress alone. 
Consequently, it is generally possible to derive a stress- 
strain relation at fixed /i in which the composition has 
been completely eliminated. This relation yields modified 
expressions for the elastic moduli with corrections c? 
in the case where Vegard's law applies. 



C. Validations 

There have been several numerical simulations of 
frame-invariant amplitude equations to date in pure 
materials^"— and alloys''^. However, quantitative val- 
idations of the results have been scarce. Here we carry 
out quantitative benchmark comparisons with known so- 
lutions to test different aspects of the method. We limit 
those comparisons to pure materials, but the conclusions 
also pertain to alloys. Validations of the alloy model, be- 
yond the derivation of the open system elastic constants 
included here, will be presented elsewhere. 



1. Asaro-Tiller-Grinfeld instability and strain relaxation 

As a first test, we study the classic Asaro-Tiller- 
Grinfeld (ATG )^^'^^ morphological instability of a uni- 
axially stressed crystal surface, which has been investi- 
gated both analytically^"— and numericallj^"— . The 
linear stability spectrum (i.e., the exponential amplifica- 
tion rate of sinusoidal perturbations of a solid-liquid in- 
terface) is known exactly for this problem and provides a 
useful quantitative basis for validation. We find that the 
dynamics defined by Eq. Q reproduces qualitatively the 
instability, but does not predict quantitatively the sta- 
bility spectrum when the wavelength is much larger than 
the interface thickness. This is because strain relaxation 
according to (HJ is diffusive, and hence artificially slow. 
If strain is not fully relaxed on the time scale that the 
instability develops, the growth or decay rate of pertur- 
bations is altered. This problem was apparently not en- 
countered in a recent PFC simulation study of the same 
instability^'^. This is possibly due to the shorter wave- 
lengths probed in this study since the PFC with diffusive 
dynamics suffers the same physical limitation. 



To overcome this limitation, we exploit the fact that 
strain relaxation in the solid can be accelerated in an am- 
plitude equation framework by making the kinetic con- 
stant r a smooth function of the density wave amplitudes. 
This allows to choose F much larger in solid than liq- 
uid. We find that, with this accelerated strain relaxation 
scheme, amplitude equations model accurately the linear 
regime of the ATG instability. We note that a differ- 
ent acceleration scheme, which has been proposed in the 
context of the PFC modeUl, consists of adding inertia to 
the equations of motion to relax the strain field propaga- 
tively, as opposed to diffusively. Although we have not 
studied this alternative here, it would be interesting to 
compare the two approaches in the future. 



2. Poly crystalline growth and grain boundary energies 

The above validation of the method for the ATG in- 
stability applies to a single crystal. For a fixed set of 
crystal axes, using the box operator or its truncation 
□j « A:^-') • V gives essentially indistinguishable results 
in the small e limit. In contrast, for a polycrystal, with 
different sets of crystal axes pointing in different direc- 
tions, the full box operator is required to make the am- 
plitude equations frame invariant, as originally proposed 
by GOS^. In order to test this frame invariance prop- 
erty, we have studied both the orientation dependence 
of solid-liquid interfaces, which controls poUycrystalline 
growth from the melt, and grain boundaries. 

As a strong test of frame invariance for polycrystalline 
growth, we have computed the solid-liquid gamma plot 
for two-dimensional hexagonal crystals, i.e. the excess 
free-energy of the solid-liquid interface, 7s;, as a func- 
tion of the angle 9 between the direction normal to the 
interface and a reference crystal axis, by two methods. 
First we keep the crystal axes fixed and vary the interface 
normal direction. Second, we rotate the crystal keeping 
the normal fixed. Both methods yield identical functions 
7s; {9) , thereby validating frame invariance. 

As a strong test of frame invariance for grain bound- 
aries, we have computed the excess free-energy of bound- 
aries, 7gh, for the whole range of misorientation for sym- 
metric tilt boundaries. These computations were car- 
ried out for both two-dimensional hexagonal and three- 
dimensional bcc ordering. In all cases, we find that 
frame invariance breaks down for large enough misori- 
entation. To illutrate this breakdown, consider symmet- 
ric tilt boundaries in two-dimensional hexagonal crystals 
with the tilt axis normal to the plane of the crystal and 
±9/2 measuring the rotation angle of each crystal from 
some fixed axis (e.g. closed packed direction), where 9 
denotes now misorientation. The initial increase of 7g{, 
with 9 is consistent with a Read-Shockley law^, as found 
previously^. However when 9 is increased further, 7g{, 
should ultimately vanish by symmetry when 9 = 60" is 
reached, since the two crystals have the same orienta- 
tion. Instead, in the amplitude equations, 7gfc continues 
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FIG. 1: Schematic representation of density waves for two 
hexagonal crystals rotated by ±0/2 where 6 is the misorienta- 
tion and the thick vertical line represents the grain boundary 
plane. The unprimed numbers label the fixed set of wavectors 
fc'-'' with J = 1 to 3 and the prime numbers label the wavec- 
tors after the rotations. Rotations are produced by spatially 
oscillating complex amplitudes. Even though the two crys- 
tals have the same orientation for the case shown here were 
9 = 60°, density waves pointing along the same direction are 
represented by different labels. 



to increase and reaches a maximum value at = GO*'. 

This unphysical feature originates from the fact that, 
even though the two crystals have the same orientation 
when 6 = 60", density waves that point along the same 
direction have different "labels" in each crystal. To see 
this, let us denote by m'^-*, u^^^ and u^'^^ the amplitudes 
of density waves with wavectors k'^^'^ — —qQ{y/?>x + y)/2, 
fc(^) = qoy, and fc^^^ — qo{^/3x — j))/2, respectively. The 
three complex conjugates u'-^^*, u^^^* and u^^^* represent 
the amplitudes of density waves pointing in opposite di- 
rections — A;(^', —k'^^\ and —k''^\ respectively, and the 
sum of all density waves represents a two-dimensional 
hexagonal crystal. For 9 — 0, both crystals are repre- 
sented by the same density waves and amplitudes such 
that 7gfc = 0. However, for 9 = GO'', the density wave 
pointing along the positive x-axis (e*"^"^) has amplitude 
u^^^ for the crystal rotated by —30", but amplitude u*^^'* 
for the crystal rotated by -|-30°, with different permuta- 
tions for the other wavectors (Fig. [TJ. Because of the ro- 
tations, those amplitudes oscillate rapidly in each crystal 
with ^ gi(-fe(3'-h9o«)-r and u^^^ ~ g-i(fe(^'-h<3oi)-? xhe 
box operator guarantees that those spatial oscillations do 
not alter the bulk properties of each crystal since the free- 
energy minimum of ([5]) is frame invariant. However, this 
frame invariance does not remove the free-energy cost 
associated with the spatially diffuse 60" rotation of each 
wavector across the grain boundary. 

This breakdown of frame invariance is analyzed in 
more detail for the simpler case of a smectic crystal in 
Appendix [b1 This analysis reveals that the grain bound- 
ary is "hidden" because the reconstructed density field 
according to Eq. ([2]) shows a perfect crystal without any 
defect. However, there is nonetheless a high free energy 
cost associated with the spatial variation of the phases 
of the complex amplitudes through the interface. This 



"hidden boundary" problem is an intrinsic limitation of 
the amplitude equations with the free-energy functional 
([5]), which extends to other crystal structures. The fact 
that grain boundaries are hidden may explain why this 
subtle issue has, to our knowledge, not been explicitly re- 
ported or analyzed in the literature so far. Extending the 
amplitude equation approach to overcome this limitation 
is an important problem for the future. 

For polycrystalline growth, this problem does not oc- 
cur because density waves do not change direction, but 
only decrease in amplitude, across the solid-liquid inter- 
face. However, a spurious grain boundary energy is cre- 
ated any time that two grains that have the same crystal 
orientation, but identical density waves labeled with dif- 
ferent amplitudes, impinge on each other. 

Despite this limitation, we note that the amplitude 
equation approach is strictly valid for low angle grain 
boundaries consisting of an array of dislocations in a 
Read-Shockley picture£2,. However, in practice, the de- 
scription remains approximately valid for misorientations 
up to about so" in the example above, or half way be- 
tween the = 0" and 9 — 60" limits where jgt vanishes 
by symmetry. Similarly, for symmetric tilt boundaries 
with a [001] tilt axis in bcc crystals, the regime of ap- 
proximate validity extends roughly up to about 45", or 
half way between = 0" and 9 = 90". Also, as demon- 
strated here, the amplitude equation approach is able to 
describe the phenomenon of "grain boundary premelt- 
ing" , which is associated with the formation of a thin 
intergranular liquid layer with a width that diverges at 
the melting point (see Refs. [23lf63l and earlier references 
therein). Since a continuous film generally only forms 
for high enough angles where dislocation cores strongly 
overlai)2^, this indicates that the amplitude equation ap- 
proach can still capture some interesting grain bound- 
ary properties beyond the Read-Shockley picture. We 
find that the dependence of the grain boundary energy 
as function of misorientation agrees quantitatively well 
with the previous PFC results^^. 



3. Grain boundary motion coupled to shear deformation 

As a last validation of the method, we have examined 
the motion of grain boundaries coupled to a shear de- 
formation. This normal motion is a generic property 
of grain boundaries that has been widely observed for 
both low and high angles^"—. It is generally faster than 
grain boundary motion associated with diffusional pro- 
cesses and hence can greatly influence the stress-driven 
evolution of polycrystalline structures. A quantitative 
understanding of this coupled motion has been obtained 
from both general theoretical considerations based on 
geometrical arguments^"— and by detailed atomistic 
simulations^"—. For low angle boundaries, this motion 
can be simply understood as the effect of Peach-Koehler 
forces on individual dislocations, which drives their mo- 
tion along a direction perpendicular to the grain bound- 
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ary plane. This motion, however, can still occur for high 
angle boundaries outside this picture. 

For perfectly coupled motion, the velocity perpendicu- 
lar to the grain boundary plane, v±, is proportional to the 
velocity v of relative translation between the two grains. 
In this case, the proportionality constant /3 = v/v± de- 
pends only on the orientations of the two crystals with 
a value determined essentially geometrically^""—. We 
find that this perfect coupled motion is well reproduced 
quantitatively by the amplitude equations for symmetric 
tilt boundaries in a two-dimensional hexagonal bicrystal 
where is analytically known. 

Atomistic simulations have shown that at high homol- 
ogous temperatures, grain boundaries can become suffi- 
ciently disordered to suppress coupled motion, which is 
superseded by a sliding motion of one grain relative to the 
other without normal motion^^. One extreme case of dis- 
order is the formation of a thin intergranular liquid film, 
associated with the aforementioned grain boundary pre- 
melting phenomeno n^^i^^ . It is clear that the formation 
of such a film will favor sliding over coupling. We show 
here that the amplitude equation approach can repro- 
duce this premelting phenomenon and the concominant 
transition from coupling to sliding when approaching the 
melting point. 



D. Outline 

The rest of this paper is organized as follows. In 
the next two sections, we construct amplitude equations 
for elemental materials in a Ginzburg-Landau spirit that 
parallels the construction of the conventional phase-field 
model. We discuss first in section II how to derive the 
"gradient-square terms" from classical DFT using a small 
gradient expansion following previous works^""— . We 
then discuss in section III how to construct the analog of 
the "double- well" potential using the equal weight ansatz 
or by deriving nonlinearities from the PFC model, which 
were obtained for bcc in Ref. \^ In section IV, still 
for elemental materials, we derive analytical expressions 
for the elastic moduli that are generally applicable to 
different crystal structures and compute their values for 
parameters of pure Fe determined previously^^ii^. In sec- 
tion V, we then extend the amplitude equations to alloys 
and derive analytical expressions for open system elastic 
constants. The various numerical validations of the re- 
sults are then presented in section VI following the same 
order in which they were summarized above. Some tech- 
nical details have been placed in appendices. 



II. GRADIENT-SQUARE TERMS 

We start with the derivation of the quadratic contribu- 
tions to the free energy density from the classical density 
functional theory of freezin g, fo llowing the procedure in 
Refs. liel - lisi (see also Refs. HiisI ). In the hquid phase 



the time- averaged particle density n{r) is spatially con- 
stant, whereas it exhibits periodic modulations in the 
solid, where the atoms have preferential positions. The 
free energy is a functional F — F[n{r)] of the atomic 
density, which can be expanded in the form 



n{r) — no + Sn{r) 
with the homogeneous density uq of the liquid and 



(9) 



*(f) 



Sn{r) 
no 



N 



(r) exp( ifc(^) • r). (10) 

3 = i 



Here, the fc^-'^ are the principal reciprocal lattice vectors; 
for bcc, these are 



[110], [101], [Oil], [110], [101], [Oil], 

[iio], [101], [Oil], [110], [101], [Oil]. 



(11) 



With the summation j = 1 . . . we write explicitly that 
we sum over the whole set of principle reciprocal lattice 
vectors, thus A^ = 12 for bcc. We point out that we 
neglect here density differences between the solid and the 
melt phase. 

The goal is to obtain an energy expression in terms of 
the density wave amplitudes ^(•'^(r). Therefore, we start 
from the known expression for the free energy change 
relative to the liquid phase, which changes due to local 
density variations. 



AF 



dfdf'Sn{r) 



no 



-C(lf-f'l) Sn{f'). 



(12) 



Here, C(r) is the direct correlation function of the liquid 
with Fourier transform 



C{q) — uq / dfC{r)exp{—ik-'r) 



with r ~ \r\^ q 
tion reads 

C{r) - 



(13) 

\k\. Therefore, the inverse transforma- 
1 



dk C{q) exp{ik ■ r). 



(27r)3no 

We also introduce the liquid structure factor 
S{q) ^ 



1 - C{q) 
First, we investigate the integral 



(14) 



(15) 



h = / dr' 



no 



C[\r^f'\) 



5n{r'). (16) 



We assume that each amplitude (density wave envelope) 
is a slowly varying function, therefore we perform a Tay- 
lor series expansion around r, 



+ - - r)kdAu^^^ (f), (17) 
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where the lower indices in the quadratic term refer to 
the vector components. This expansion is inserted into 
the above integral expression. The contribution from the 
(5-function gives readily 

5{r~r') 1 

dr'-^ '-5n{r') = —Snlr). (18) 

no no 

Second, the contribution from the r' independent term 
in the series expansion (fT7| leads to the integral expres- 
sion 



r — r X 



N 

-no^u(^)(r)exp(ifc(^) -r) / dr'C{\ 
X exp[ifc'-'-' • (r" — r)] 

N 

= - ^ w(^) (r) exp(ifc(^') • r)C{q), 



since C{q) — C{q)* due to the inversion invariance of 
C(|f|) (the star denotes complex conjugation). 

Next, we note that the linear term in Eq. PT)) does 
not contribute to the remaining part of /i, since we as- 
sume that all k-vectors are at the (highest) peak of the 
structure factor, i.e. C'{qo) = 0; in fact, this term is 



^no 



^ exp(ifc(^') • r) J dr'{r' - r) ■ [Vu^^'^ (f)] x 

xC(|r-f'|)exp[ifc(j') • {r' -f)]. 



We therefore inspect the integral (with the notation q' 
\k'\) 

no / c?r '(r ' — r)C(|r — r 'I) exp[ifc • (r ' — r)] 



no j dr'r'C{\r'\)eyi^[ik-r'] 
1 



(27r) 



dr'r' I dk'C{q')ejip[i{k + k')-r' 



/ df' / dk'Ciq')exp[i{k + k') 



(27r)3 ^ 

= -i'^kj dk'C{q')5{k + k') 

= -iVj:C{qo)={). 

Finally, we look at the term that arises from the 
quadratic term in the Taylor expansion. First, we show 

no I dfC(|r=1) exp(^fc • r)nr, = ~C"{qo)kik^ (19) 



where kj = kj / qo are the normalized components of the 
vector fc; the integral is performed on the entire space. 
Thus we get 



no 



J dfCQ 



f\) exp(ik ■ r)firj 



(27r)3 ' / *'^('?')'2xp[i(fc + fc') • 
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(27r)3 dkidkj 



df / dk'C{q')exp[i{k + k') -f] 



dkidkj 
dkidkj 



-C'iqo)'^ 
% 



dk'C{q')5{k + k') 
C{qo) 

-C"{qo)kik 



since in the last steps C'{qo) = 0. With these prerequi- 
sites, the last remaining term in Ii becomes 



1 ^ 

-no^exp(ifc(-') • r) {dkdiu''^\f) 



X / df'{f' - f)k{f' - r)iC{\r' - r\) exp[ifc(^') • (f - r)] 



1 ^ 

= - ^ exp(ifc(^) • f) (dkdiu^'\f)) C"(go)A;^'fcp' 
i=i 

1 ^ 

= - ^ exp(ifc(^) • r)C"(go)(fc*^'^ • '^fu'^^\f) 

Altogether, we get 
w -, 

h = V — — u^^Vr) exp(iA:(j) • r) (20) 

+ -^exp(ifc(j') • f}C"{qo){k'^'^ ■ V)2wW(r), 

thus the expression for the free energy becomes 
keT 



df5n{f)Ii 



(21) 



Now we assume a separation of scales, which means that 
the scale over which the amplitudes vary is much longer 
than the atomic spacing ^/qo- Hence we get for a 
"slow" function gg 

dr exp(z(fc(^) + fc«) • f)g,{f) J df <5,„,,^,,„ ^^^.(r). 

(22) 



Therefore 



AF = 



noksT 
2 

N 



+ -^C"{qo)u^^^*{f){U^'^ • V)2u(j')(f) 



(23) 



where we used that it^'^ = u^^^* if fc''-' = —k^^\ since 
the density is real. We can also integrate the last term 
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by part, assuming that the boundary terms do not con- 
tribute for appropriate boundary conditions and obtain 
finally 



AF 



2 

N 



df 



1 ^ 



Siqo) 



uO-)u(i)* 



|EC"(go)|(fc(^')-V> 



.0)1 2 



J = l 



(24) 



As pointed out in Refs. ISOHSfl the operator fc'^^ • V 
violates the rotational invariance of the functional. As 
already mentioned, the proper renormalization is to re- 
place it by the "box operator" 

The properties of this operator and the consequences for 
the model will be discussed in detail in section IVIBI 

The derivation of the free energy allows to link the 
phenomenological parameter ^ in Eq. ([5]) to physical pa- 
rameters via the relation C — —n^kBTC" [qo) / A. 



III. DOUBLE WELL POTENTIAL 

We now discuss the derivation of higher order cubic 
and quartic nonlinear terms. With the addition of those 
nonlinearities, the bulk free-energy density has the form 
of the standard "double- well" potential of the phase-field 
model. As already mentioned in Section HI nonlineari- 
ties can be obtained in a Ginzburg-landau spirit using 
the equal weight ansatz of Shih et ali^, or derived from 
the PFC modelS^i^i^. Here we give the results of both 
methods for bcc and 2-d hexagonal cases. 

Due to the assumption of a scale separation, orthog- 
onality demands that only higher nonlinearities which 
form a closed polygon of reciprocal lattice vectors con- 
tribute. The addition of cubic and quartic terms gives 
therefore the general expression 



AF = / dr 



Silo) ^ 



,0)1 2 



(26) 



C"(go) 



N 



,0)1 2 



ijk 



+a4 E '^Ofc'"^'^^^''^"^'''^"^'^'^o,fc(*)+feO)+fe('»)+fe(') 

ijkl 



In the summation over three or four wave vectors in the 
cubic and quartic terms the summation is normalized 
such that permuting indices that correspond to equiva- 
lent sets of fc- vectors are only counted once. Using an 



equal weight ansatz, one obtains^ 



and 



- 1/ 



0.3 



= 1/27 



a4 



24 
S{qo)us' 
12 

5(qo)ur 



(27) 

(28) 
(29) 



where is the amplitude of all density waves in the solid 
phase^i^. The connection to a conventional phase field 
model with a double well potential becomes obvious in 
the isotropic approximation, where all density wave am- 
plitudes are equal, u'^' = u, and assumed to be real. 
Then the local part of the free energy density ((26)) be- 
comes 



noksT 12 
2 S{qo)u': 



(30) 



which obviously has energetically equivalent minima for 
the bulk states w = and u = u^. 

Alternatively, the coefficients of the higher order non- 
linearities can be derived using a multiscale analysis of 
a phase field crystal model, yielding different expressions 
for the quartic coefficients only^: 



^ijkl ~ 



only two wave vectors fc(^) and -fc(^) 



all other quartic 



(31) 

So far, the minima of the free energy correspond to 
solid and liquid phases which are energetically equiva- 
lent. A deviation from the melting temperature Tm fa- 
vors one or the other phase, and this is achieved here by 
introducing a tilt term of the form 



T, 



M 



(32) 



with the latent heat L and a coupling function (anal- 
ogous to the standard "phase field") that has value 1 in 
the solid and in the liquid; the function ip should be sta- 
tionary there in order not to shift the bulk states w*^^^ = 
and u^^^ = Us- A particular choice is 



0({u(^)}) 



1 



N 

^EmM^'^iV"^ 



with 



h{x) = x^(3 - 2x) 



h{x) = x{i - 2V^). 



(33) 



(34) 



(35) 



For the simulations shown here we used the first choice, 
Eq. Other choices are possible as in the conven- 

tional phase-field approach. 
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For the purpose of numerical implementation, it is use- 
ful to rewrite the free energy in a dimensionless version. 
We therefore introduce a new small parameter 



24 



Siqo)C"{qo)ql' 
a dimensionless amplitude 

and a dimensionless (slow) lengthscale 



(36) 



(37) 



(38) 



We note that e is defined differently here for bcc than 
in Ref. to eliminate all PFC parameters from the 
amplitude equations. Then the free energy becomes for 
the phase field crystal model 



\-j=i j=i 
r (Nil \ ^ NI2 

+2^iio^iio"^ioi^ioi + 2^110^110^101^101 
+2^110^011^011^110 + 2^tlo^oii^oii^iio 

-g l^oii^ioi^Iio + ^ou^toi^iio 

+^011^110^^01 + ^011^110^101 
+"^011^110^101 + ^011^110^101 



+^011^101^110 + ^Oll^loI^llO 



with 



Fn = - 



C"iq„)q^'uis 



-1 ,2^-1/2 



(39) 



(40) 



Here, the box operator is defined as 



le 



1/2 



% 



le-i/^D,, (41) 



where the nabla operator V acts on the variable R. 
The amplitudes A*^^^ are defined as functions of the (di- 
mensionless) "slow" scale X, which is much larger than 
the atomic spacing for e <C 1. The summation N/2 
in the expression l\'39\i expresses that we sum only over 
the contributions from independent density waves, since 
^0)* = ^0) , jt means e.g. for bcc that we sum only over 
the six principal reciprocal lattice vectors [110], [101], 
[Oil], [110], [101], [Oil]. 



Similarly, the thermal tilt ((32)) becomes 



Tm 



with 



N/2 



(42) 



(43) 



j=i 



The reconstructed density becomes 
n{R) ~ TiQ 



l + ..f:A(^)(i?)exp^^^^'^-^^ 



i=i 



no 



N 



1 + u^^^ (R) exp 
i=i 



£l/2 

ifc(j) • r\ 

£1/2 



(44) 



The dynamical evolution of the nonconserved order pa- 
rameters is described by 



dt 



SF 



(45) 



As long as we are interested only in equilibriums proper- 
ties, the choice Kj > (or correspondingly F in Eq. (|3])) 
is not important. The formulation of the dynamics will 
be discussed in detail in Section IVI Al 

Similarly, for a two-dimensional hexagonal model we 
get the (rescaled) free energy from a phase field crystal 
model (see Appendix |X] for details) 

(N/2 2 1 

I ]=i j=i 



-^^(itiSi^ + iiisia) 

/7V/2 \ ^ N/2 

\J=1 / J=l 



(46) 



where we have N = 6 principal reciprocal lattice vectors, 
see Eq. (jASp . Assuming an equal weight ansatz, we again 
get slightly different quartic terms: 

f N/2 2 1 

F^D - F^D / dR\j2\d,A,\ +gE^^^.* 
I j=i j=i 

(n/2 \ ^ N/2 

+^ EM* 

\j=i / j=i 



(47) 



IV. ELASTICITY 

The amplitude equations naturally incorporate linear 
elasticity. Here we derive an analytical expressions for 
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the elastic constants that are generally applicable to dif- 
ferent crystal structures by summing the contributions of 
different sets of crystal density waves. Those expressions 
reduce to linear elasticity for two-dimensional hexagonal 
crystals and yield reasonable predictions of the elastic 
moduli for bcc. Both hexagons and bcc can be repre- 
sented by only the principal set of reciprocal lattice vec- 
tors. In contrast, for the PFC model of fee structuresSS, 
the present analysis shows that the addition of a second 
set of density waves is required to obtain a non- vanishing 
tetragonal shear modulus. Here we restrict our attention 
to linear behavior elastic behavior. Nonlinear elasticity 
has also been shown to be analytically treatable in an 
amplitude equation framework^. 
We start from a density field 



dno{r} 



no 



r) exp (ifc(j') • f) (48) 



of a stress free system. The amplitude u^^'^-' (r) may de- 
scribe a pure solid or a system consisting both of solid 
and liquid parts, and we only assume that the amplitude 
varies slowly in space. If e.g. an external stress is applied 
to the system, the atoms are displaced and take new po- 
sitions. In other words, an atom that was previously 
located at r is now at r + u{r) , with u being the displace- 
ment field. Ignoring density changes due to the strain, 
this leads to the condition 5n[f + u{r)) = 5no{r)- Since 
we assumed that the amplitudes themselves are slowly 

(J'°)(r). 



u 



varying functions, we can assume u^^''^^ {r + u) 
Then we obtain for the density field of the stressed sam- 
ple 

5n{r) — uq u^-'''^\f) exp(— ifc'-'-' • u{r)) exp(ifc^^' • f). 

(49) 

This means that the material is now described by new 
density wave amplitudes 

uO)(f) = M(j''°)(f) exp(-ifc(^') • u{r)), (50) 

which are oscillating functions (in a solid) , and the wave- 
length of the oscillation is shorter for higher strains. 

The principal reciprocal lattice vectors obey the fol- 
lowing orthogonality relation-^ 



N 



(51) 



where N is the number of principal reciprocal lattice vec- 
tors, d the spatial dimension and qo = which is 
equal for all principal reciprocal lattice vectors. Here, 
the lower index denotes the component of a vector. It 
is straightforward to check that this relation holds for 
the two present cases of interest, the three-dimensional 
bcc lattice, Eq. ([TT|) and the two-dimensional hexagonal 
lattice, Eq. ^M^. 



We can therefore extract the displacement field from 
the amplitudes by taking the (complex) logarithm 



Ln 



M(j)(f) 

\uU){r)\ 



= ■ u{r) + 2TTim{r) + icj)^. (52) 



Here, we used the symbolic notation of an integer "wind- 
ing number" m(r) which stems from the fact that the 
imaginary part of the complex logarithm. 



Ln z =^ In Izl 



I arg z, 



(53) 



is defined only up to arbitrary shifts by 27ri, which can 
be absorbed in the complex phase argz. Since we are 
mainly interested in derivatives of the above expression, 
this additional integer term usually disappears anyway. 
It becomes only relevant if defects are present in the sys- 
tem. The last term (/iq is the phase of the slow ampli- 
tude u^^'^\ which may result from a rigid body transla- 
tion, and which is therefore constant. Another possibility 
would be a rigid body rotation, and then (/)o is not con- 
stant. However, this term does not contribute to the final 
result and we therefore ignore it here. 

Differentiation and summation over all fields gives, to- 
gether with Eq. 



N 



N 



J^fcp^a^Ln^^ = -z(a„u„(f))5:fcPfcW 



'i—qoOmUi[r). 



We therefore arrive at the following expression for the 
(infinitesimal) strain tensor e/m — {diUm + dmUi)/2: 



i d 



iw(f)l' 



(54) 



By straightforward algebraic manipulations we obtain 
the alternative representation 




L(fcWa™-ffc2^50uW I , (55) 



where 3(-) denotes the imaginary part. 

In the expression for the free energy, only the gradient 
term changes if a solid phase is displaced. The local terms 
remain unchanged, since by construction they consist of 
density waves amplitude products with closed polygons 
of principal reciprocal lattice vectors. Therefore, in any 
product like w^^^m^^^u*^^^ (for the hexagonal system, see 
Appendix \K\ for details), we have 



^(1,0)^(2.0)^(3,0) ^ 

X exp[-i(fc(i) -f 

^(1,0)^(2,0)^(3,0) 



and it is therefore sufficient to inspect the gradient terms. 
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The "kinetic" part of the free energy is 

N 



1 noksT 



C"{qo) J dr^KfcW-Vh 



a)i2 



(56) 



Notice that we dropped here the correction term from 
the box operator. The reason is that this term gives only 
a higher order correction, and for elastic deformations 
(which are always assumed to be long- wave distortions), 
the additional term is negligible. Let us assume for sim- 
plicity that the original amplitudes u^^'^^^ are real. Then 
we get 

|(fc(^).V)M«P = (^P(9m(^'°))) V(fcp^/fc^)u(^'"^9mm)' . 

(57) 

The first term is the usual interfacial energy (the same 
as for the undeformed state), and the second term, which 
is quadratic in the distortions, the elastic energy. Notice 
that we do not get a term that is linear in the strain, 
and therefore we do not have a surface stress term in the 
model. 

We get explicitly for the elastic term 



InakBTC'iqo) 



N 



(58) 



We can therefore identify the elastic constants 

N 



"0^^^ ^''(go) V A-^^') k^^^ k'^^KiU.oy 



(59) 

which are obviously invariant under pairwise permuta- 
tions of indices. Here we used 



Cijkl^kl 



(60) 



and the usual expression for the elastic free energy den- 
sity. 



/e, 



(61) 



Notice that the expression for the elastic constants au- 
tomatically reflects the correct crystallographic symme- 
tries. 

First, we note that in the "liquid" phase all elastic 
constants become zero; this is a consequence of the fact 
that we skipped contributions from density changes. 

For the bcc case with cubic symmetry we obtain 



2 \ii=j^k = l 

Cijki = C*" (90)90^5 ^ { 1 two distinct index pairs 

else 

(62) 

We can calculate explicitly the predicted values for the 
elastic constants for bcc iron at the melting point. Using 
the parameters given in Ref. which are summarized 
in table HI we obtain (in Voigt notation) 



Cii 
C12 



C22 
C23 



C33 = 90 GPa 
C44 = 45 GPa. 



(63) 
(64) 



Given that the theory only includes one set of density 
waves, those values are reasonably good (see Ref. [2^ for 
a quantitative comparison with MD results) . 

The two-mode PFC model for fee structures couples 
two different sets of crystal density waves corresponding 
to (111) and (200) reciprocal lattice vectors25. Using the 
analytical expression for the elastic constants, Eq. ([5^ . 
it is straightforward to work out that the (111) set yields 
equal elastic moduli C12 = C'22 = C44, and hence van- 
ishing tetragonal modulus (Cn — Ci2)/2. However this 
unphysical feature is cured by the addition of the (200) 
set that brings an additional finite contribution to Cn 
and vanishing contributions to C12 and C44. Therefore, 
with both the (111) and (200) sets present, Cu > C12 
and C12 = C44, such that the tetragonal shear modulus 
is finite as physically desired. Furthermore, the expres- 
sions for the elastic moduli predicted by Eq. ((59)) agree 
with those derived in Ref. |25j by a brute force calcu- 
lation of quadratic contributions to the PFC two-mode 
free-energy functional for different lattice distortions. 

Next, from the equilibrium conditions we know that 
the free energy F has to be minimized with respect to all 
degrees of freedom. For fixed interface position this im- 
plies that in particular F has to minimized with respect 
to the elastic displacements, which appear only in the 
elastic contribution fei- Therefore, we get 6F/Sui — 
and consequently 



dxi 



(65) 



which are the usual elastic equations. 

Finally, we note that the model allows for deformations 
that are not contained in the standard theory of linear 
elasticity. The reason is that the displacement vector has 
only d components, but we have have N/2 independent 
amplitudes. This means that not all possible amplitudes 
can be represented in the form ([50]) with real amplitudes 
yO'.o) Qf^]-^g "undeformed" crystal. The remaining iV/2 — 
d degrees of freedom correspond to atomic "shuffles" i.e. 
rearrangements within each unit cell. 



V. ALLOYS AND VEGARD'S LAW 

We now consider the extension of the amplitude equa- 
tions to binary alloys. As discussed in section I, the 
present extension has the advantage that it can interpo- 
late between the thermodynamic properties of the solid 
and liquid phases, which can be described in principle by 
arbitrary free-energy/composition curves. For concrete- 
ness, we consider here a binary alloy in the dilute regime. 

The impurities are introduced variationally using a 
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TABLE I: Parameters for bcc iron, see Ref. |43 for details, as obtained from MD simulations^^. 



go 


S(go) 


C"{qo) Us no 


Tm 


L 


2.985 ■ 10"m" 


3.01 


-10.4 ■ lO^^^Ti^ 0.72 0.0765 ■ lO^V"^ 


1773 K 1.968 


10^ Jm"^ 



new free energy term 
RTm 



Fc+Fj 



df 



c In c — c) + (j)Ae c + L 



Tm 



(66) 

which includes the previous temperature coupling and 
the phase field (p as defined in Eq. It contains the 

molar volume vq and the ideal gas constant R. 

The evolution equation for the amplitudes is the same 
as before, 



SF 



dt 5u(3> ' 

where we get now the additional term 



(67) 



(i) 



The diffusion equation follows from 



dt 



D 



^^F 
cV — 

RTm oc 



Nul' 
(68) 



(69) 



with a diffusion coefficient that can be different in solid 
and liquid, and gives therefore 



with 



dc 

— = V • [DVc - DbcVS 
ot 



«oAe 
RTm 



(70) 



(71) 



These equations describe in equilibrium a phase dia- 
gram with straight solidus and liquidus lines. The parti- 
tion coefficient is given by 



k = exp(&), 
and the liquidus slope is 

VqL 



(72) 



(73) 



Additionally, we can take into account that the equi- 
librium lattice constant changes with impurity concen- 
tration. For a linear dependence (Vegard's law^i^) we 
can change the box operator to the gauge invariant form, 
as discussed in section |T] 



U, ^ 



(74) 



where a is proportional to the expansion coefficient. As- 
suming a long- wave modulation of the concentration, we 



can ignore the higher order correction term in the box 
operator. 

For simplicity, consider a deformed solid with 



,0) 



Us exp 



[-ik^3) -uir)). 



(75) 



Then we already know that all the local terms in the free 
energy functional remain invariant under elastic defor- 
mations, thus it is sufficient to look at the gradient term. 
Thus we get for each wave vector 



iac\u 



0) = -iu^i\qaUj;^k\''^diu,^ - «c) (76) 



0) 



Thus the expression in the free energy functional, |(fc' 

V -f lac)?/^-'-' P, is minimized for acq^^ — k!'mk\^^diUm- 

Symmetrization gives immediately acq^ — km ki eim- 
Using the orthogonality theorem ([5T|l we therefore get 
the relative lattice expansion 



ta = dacqQ 



(77) 



One can readily check (e.g. for the hexagonal lattice 
(|A5p ) that all diagonal elements of the strain tensor are 
equal and that the off-diagonal elements vanish, thus we 
get a dilatational stress free eigenstrain 



and the elastic energy density becomes 

fel = -j^CijkliUj - ^°j){^ki - ^li) 



(78) 



(79) 



with the same elastic constants as before. 

Notice that the lattice dilatation via the modification 
of the box operator Eq. ([71]) also affects the impurity 
diffusion, since it introduces another concentration de- 
pendent term in the free energy. The "kinetic" part of 
the free energy 



AFi.. 



mkBTC"{qo 



N/2 

■E 



dr 



(□j -f iac) u 



gives a contribution to the chemical potential, 
5Fk _ nokBTC'iqo) 



(80) 



Sc 2 

N/2 

X ^ [ - 2ak^^'> ■ 9(u(J')m(J>) - aqo^n{u^^'>V'^u(^^'>*) 
+2a^cu^^'>u'^^'>*] . (81) 
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For a deformed solid, Eq. ((75|) . we obtain 

n^kBTC"{q^) r , N 
IJ-k = 7; [ - au^qo — diui 



1 



N 



+ -auiqo — i^m)i^^Ul) + Na'cuil (82) 



where we used the orthogonahty relation (|5ip . As before, 
the second term stems from the higher order term in the 
box operator, and as long as the distortions are small, 
diUi <C 1, it can be neglected (in fact, it corresponds 
to the nonlinear contribution in the full strain tensor, 
fifc = {dkUi + diUk + diUidkUi)/2, which is relevant for 
its rotational invariance). This is typically the case if the 
eigenstrain is small, ucqq^ 1. For a deviation from the 
equilibrium strain, e-ij = e^^ 



Seij, wc therefore obtain 



fJ-k 



nokBTC'iqo) 



N 



2 -au;qo — Seu 
= -go Vc,jfefe(e,y - e°-) 



(83) 



where we used the general expressions (l59l) and (|60p . No- 
tice that the same expression can also be derived directly 
from Eq. ([7^ and fik — dfei/dc. 

If the material is stretched. Sen > or akk > 0, the 
chemical potential is reduced, /Xfc < for a positive Veg- 
ard coefficient a > (the material extends if it contains 
impurities). Thus a flow j —cV^k sets in to locally in- 
crease the concentration. Hence, the impurity concentra- 
tion is affected by the local volume change. The diffusion 
equation becomes therefore instead of Eq. ([70)) 



dc 
dt 



DWc - DbcWc/) + D 



RTm 



(84) 



with the expression for /x^ being given by Eq. (jSTj). 

It is straightforward to generalize the model to cases 
with more complicated phase diagrams, eventually also 
using thermodynamic databases. Similarly, the extension 
towards nonlinear or nonisotropic lattice expansions and 
the inclusion of thermal expansion using a temperature 
dependent expansion coefficient is straightforward. 



A. Open-system elastic constants 

The elastic constants Cijki describe the material stiff- 
ness for fixed composition. As mentioned above, stretch- 
ing of the material leads to an increased solute concen- 
tration and induces elastic relaxation. Therefore, for fast 
diffusing elements and slow deformations, the material 
seems to be effectively softer, because it compensates 
the elastic deformation by e.g. filling the interstitial po- 
sitions with impurities (increase of the stress-free strain 
e°j). Formally, this leads for fixed chemical potential 
to the definition of the open system elastic constants^, 
which can be calculated readily using the expressions 
given above. 



The chemical potential becomes in the solid with the 
small strain approximation, see Eqs. ([M]) and 



SF 
1^ 



RT, 



M 
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In c + Ae - 



9o '^'^kk , 



(85) 



and therefore the local concentration as function of stress 



Co exp 



^090 ^ 
RTm 



OLGkk 
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RTm 



aakk 



with 



cq = exp 



vo% ^ 
RTm 



(/i-Ae) 



(86) 



(87) 



For positive Vegard coefficient, a > 0, Eq. (|86|) expresses 
the local concentration increase under tension, a^k > 0. 
Therefore, the stress-strain relation aij — CijkiX^ki — £/!;) 
becomes implicit through the stress free strain (l78l) 



^ij ^ijkl 



- aco^o "fe' 



RTMql 



Op'SklCTmrn 



This relation can be inverted, 

o-y = c*jki[<^ki -eH(co)], 



(89) 



which defines the open-system elastic constants cl-^i- 
For a general case of cubic symmetry, with only 



--xxxx-i ^xyxy 



and 



-■xxyy 



being independent elastic con- 



stants, we obtain from Eqs. (IH5)) and 

* ^xxxx ^~ '^X.V i^xxxx ~l~ ^^xxyy^i.^: 



yy J \^xxxx 



^xxyy) 



1 ~l~ i^xxxx H~ '^(^xxyy^ 

^xxyy ~t~ \Tj {Cxxxx ~\~ ^^xxyy^ij^xxyy 



^xyxy 



1 ^~ ^W'^i.^xxxx H" '^Cxxyy^ 

^xyxy 

(90) 

with X s-iid 77 being defined as in Ref. [52I 

Co^'O 



V = c,qo^- (92) 

Interestingly, the combination / :— Cxxxx — Cxxyy — '^Cxyxy 
remains the same for the open system, i.e. /* = /. No- 
tice that / = is the condition for isotropy of a cu- 
bic systenJi^ and since the impurities change the lat- 
tice constant uniformly in all directions, the solute does 
not destroy the isotropy. In that case, the material can 
be characterized through two elastic constants, e.g. the 
shear modulus and the bulk modulus K. From the 
general expressions ([50)) we obtain for / = the isotropic 
open system constants 



1 1 



K* 



K 



(93) 
(94) 
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These expressions match exactly the prediction in 
Ref. 152]. The first equation reflects that a pure shear 
does not change the volume, and therefore no concentra- 
tion change occurs. As expected, the open system bulk 
modulus is smaller than for fixed concentration. We note 
that this reduction occurs also for a < 0, since the cor- 
rection is quadratic in a. 



VI. VALIDATIONS 

A. The Asaro-Tiller-Grinfeld instability 

The Asaro-Tiller-Grinfeld (ATG) instability is a 
morphological instability of a uniaxially strained 
surfac o^'^'^^" — . The development of corrugations due to 
a reshuffling of material reduces the total energy for long- 
wave perturbations. Here, local melting and solidification 
at the interface leads to the development of the insta- 
bility. In two dimensions, the shape of the interface is 
described by the profile (see Fig. [2]) 



x{y) — A sin ky 



(95) 



and subjected to a tensile or compressive stress ctq along 
the interface. Since we assume that the melt phase is 
stress free, the normal and shear stresses vanish in the 
solid at the interface. Then the chemical potential dif- 
ference at the interface between the solid and the melt 
becomes in a sharp interface picture^ 



(96) 



with the atomic volume 51 and the interface curvature 
K. Since the the solid- melt interfacial energy j{9) is 
anisotropic, the stiffness jc« = {9) appears here. 

It triggers interface evolution via a melting- 
solidification process, and the interface normal velocity 
is given by 



M 



A/i 



(97) 



with the kinetic coefficient M of the interface kinetics. 

The evolution of the interface leads to a time- 
dependend amplitude (in the framework of a linear sta- 
bility analysis) A = Aq exp wt. A sharp interface calcula- 
tion predicts for isotropic elasticity in a two-dimensional 
plane-strain situation the spectrum 



M 



- e-f,sie = 0) 



(98) 



where M is the kinetic coefficient of the melting and 
solidification process, and the (planar) interface normal 
direction is assumed to correspond to 6* = 0. E and 
!/ are Young's modulus and Poisson ratio respectively. 
Here, the first term accounts for the elastic destabiliza- 
tion, whereas the second term describes the stabilization 



due surface energy. The above spectrum sets a charac- 
teristic lengthscale, the Grinfeld length. 



Lg = 



2(1-^.2)^2 



(99) 



We model this process using the two-dimensional am- 
plitude equations for a hexagonal system and also use 
this notation (see Appendix E|) . We note that the six- 
fold symmetry induces elastic isotropy, whereas the in- 
terfacial energy remains anisotropic (see also Fig. [S] in 
Section [VI B I) . For the amplitude equations we use relax- 
ation equations of the type 



dt 



-K 



6F 



(100) 



with kinetic coefficients Kj. These equations do not lead 
to the same sharp interface limit as used for the ATG 
spectrum above if all Kj = K are constants. The reason 
is that the motion of the interface occurs on the same 
timescale as the relaxation of the elastic degrees of free- 
dom. This is problematic especially in the long wave 
limit, because the range of the elastic distortion is the 
same as the wavelength of the interface corrugation, and 
the elastic fields have to adjust via a diffusive process. 
In reality, however, the interface motion is slow in com- 
parison to the sound speed (which sets the true scale for 
the elastic relaxation), and therefore the assumption of 
static elasticity as in Eq. ([98)) is appropriate. 

To obtain the same behavior with the amplitude equa- 
tions, we use a kinetic coefficient that depends on the 
amplitudes: In the solid, the kinetic coefficient Kj — Kg 
is high and low in the liquid, Kj = Ki. In between, the 
coefficients are interpolated. 



Kj = hjKs + (1 - hj)Ki. 



(101) 



with interpolation functions hj that have value 1 in the 
solid and in the liquid. This is the aforementioned de- 
pendence of the kinetic coefficient F (in the notation of 
Section|l| on the local values of the amplitudes. For sharp 
interfaces, i.e. the width of the diffuse interfaces ~ e 



-1/2 
2D 



being small in comparison to the wavelength of the per- 
turbation 1/fc, the precise choice of the interpolation 
is not crucial. In particular, we used hj = ft,(|A°/^°p), 
with h{x) being given by Eq. (j34p . Since the motion of 
the interface is basically determined by the smaller of the 
coefficients K^^ Ki, we can therefore get a slow motion of 
the interface, whereas the elastic relaxation in the solid is 
sufficiently fast, since it is determined by K^. In the limit 
Ki/Kg we therefore recover the case of quasistatic 
elasticity. 

Since the kinetic coefficient M cannot easily be ex- 
pressed in terms of the mobilities Kg and Ki, we first 
investigated the decay of capillary waves without elas- 
tic effects, i.e. without application of an external stress. 
From the decay rate and the spectrum (|98p we therefore 
extract the value of M for given values of Kg, Ki. Next, 
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we apply additionally an elastic deformation tangentially 
to the interface and measure the modified amplitude evo- 
lution. Snapshots of the temporal evolution in the unsta- 
ble regime are shown in [5J The simulations are started 
with a small initial amplitude /cAq ~ 0.11, which grows 
here for kLc — 0.84. The interface thickness ^ used in 
the simulations is related to the wavelength of the per- 
turbation by fc^ ~ 0.3. 

We use a straightforward real space discretization of 
the amplitude equations, and fixed boundary conditions 
in the direction perpendicular to the interface {x direc- 
tion, see Fig. [5]): In the liquid (x = 0) the amplitudes are 
fixed to zero, whereas at the right interface {x — X) we 
have 

Ai{X,y) = A.expHfc^'^ -"0(^,2/)], 
A2(X,y) = -A,expHA:(2).^o(X,y)], 
A:i(X,y) = A,exp[-ifc(3) .^^(^^y)]^ 

with the homogeneous displacement field for the planar 
front 



xel, 

y4y 



which depends on the homogeneous strains and Eyy. 
Notice that for a homogeneous nonhydrostatic stress (Tq 
the system needs to be strained in both directions ac- 
cording to 



E 



yy 



E 



-(Jo- 



in the other direction along the interface (y direction), 
we use quasiperiodic boundary conditions: 

Aj{x,Y) = Aj{x,0)ey.^[-iU'^^ ■ Auq] 

where Y is the system size in this direction, and the dis- 
placement jump Alio = U{){x,Y) — Mo(x, 0). The wave- 
length of the perturbation therefore has to fit into this 
periodic interval. Instead of changing the wavelength to 
scan the spectrum of the ATG instability, we vary the 
stress tJo; this also has the advantage that the scale sep- 
aration between the different geometrical length scales 
does not change. 

Due to the energy increase of the solid through the 
mechanical load, it is no longer in equilibrium with the 
melt at T = Tm- It is therefore convenient to suppress 
the planar front motion by a slight undercooling, since 
we focus here on the development of the instability. 

For the calculation of the spectrum, we use only the 
first linear regime of the amplitude evolution (after an 
initial stage where the interfaces adjust to the proper 
profiles). The results are shown in Fig. [31 together with 
a comparison to the sharp interface prediction (|98|. We 
clearly see that in the limit Ki/Kg the spectrum 
agrees well with the analytical theory. In particular, for 
kLc w 1 we get w ss 0; then the interface motion is slow. 



FIG. 2: Interface evolution of the ATG instability for kLc = 
0.84 with Ki/Ks = 1 and Kb^2d = 1. The "atomic spacing" is 
determined through the parameter e2D and is relevant only for 
the reconstruction of the density waves but does not influence 
the "macroscopic" evolution of the interface. Here e2D = 
0.5 and periodic boundary conditions are used in the vertical 
direction. The time in the snapshots is from top to bottom 
0, 12000, and 20000, respectively. At a nonhnear stage of the 
instability, deep grooves form in the solid to reduce the elastic 
energy, which can lead to fracture -4i22. 
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FIG. 3: Spectrum of the ATG instability. The sohd hne is the 
analytical sharp interface prediction Eq. (|98|l . the dots show 
the numerical results obtained from the amplitude equations, 
for two different mobility ratios. For small Ki/Ks we recover 
the sharp interface limit because the strain relaxes fast in 
comparison to the interface motion. 








and therefore even for equal mobility in solid and melt 
the elastic fields can adjust fast enough. Hence curves 
for different mobility ratios intersect all at kLc = 1. 

Since the rotation of the grains is not important here, 
the higher order correction in the box operator can be ne- 
glected. We checked numerically that for €20 ~ 0.1 the 
correction term gives only negligible modifications of the 
results. Then the small parameter €20 appears in the 
free energy functional Eq. (|A22p and in the amplitude 
equations (jlOOp only as multiplicative constant and can 
be absorbed in the kinetic coefficients, thus the descrip- 
tion is entirely on the slow scale R. The parameter €20 
comes in only via the reconstruction of the density waves 
according to Eq. (jA18|) . In particular, for a solid that 
is not rotated the amplitudes are constant in the solid 
(without strain) or vary only gently if a strain is present, 
and therefore the amplitudes can be discretized on a scale 
that is independent of the "atomic" resolution. In this 
sense, the computational efficiency of the model is not 
inferior to a conventional phase field model, apart from 
the fact that more than one parameter is needed. On 
the other hand, the description automatically contains 
elasticity, which would require a separate treatment in a 
conventional model, see e. g. Ref. nt 



B. Crystal-melt interfacial free energies and 
polycrystalline growth 



A central part of the theory is the box operator, which 
generalizes the gradient term fc^^' • V of a more conven- 
tional Ginzburg-Landau theory to a rotational invariant 



FIG. 4: Schematic illustration of different rotations of a crys- 
tal and/or its interface contour, (a) Original state with ampli- 
tude A{R). (b) The contour of the crystal is the same, but the 
lattice is rotated; the amplitude is A{R) exp{ik^'M.R/e^^^); 
(c) The contour is the same, but the lattice is rotated in the 
opposite direction; ^(_R) exp(ife^M^_R/£^''^). (d) The contour 
is rotated, but the lattice orientation is the same as in (a); 
the amplitude is yl(R_R). (e) Both the contour and the lattice 
are rotated; therefore this state is equivalent to (a) and the 
amplitude is ACRE.) exp{ik'^MR/e^^'^). 



form. 



1/2 



(102) 



which introduces a higher order correction. We use here 
the dimensionless representation introduced in Eqs. ([36])- 
(UH). Since it also brings higher order derivatives, a nu- 
merical treatment becomes computationally more costly 
in an explicit scheme, since then the relaxation timesteps 
have to be rather short. 

In the following, we pick a particular k-vector, and 
drop therefore the subscript j. In a pure solid phase 
with a spatially constant amplitude A, we have DA = 0. 
However, this relation also holds if we describe a solid in 
a rotated state. Namely, if A{R) is an amplitude field, 
a crystal with the same shape, but rotated grain struc- 
ture, as sketched in the transition from Fig. to b, is 
described by 



A+{R) = A{R)exp 



(103) 



where the dagger f denotes transposition. Here, M = 
R — I, where I is the unity matrix and R an orthogonal 
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rotation matrix. In two dimensions, R has therefore the 
structure 



R 



cos 9 sin 9 



— sm ( 



(104) 



grain with a lattice that is rotated in opposite direction 
(Fig. dJ:) . To obtain the relation (|108p we first note that 
V^(R-R) = R^Vyi|j^j^. Furthermore, the Laplace oper- 
ator is rotational invariant, i.e. V^^(R^) ~ 
Then we get 



The reason is that the rotated lattice structure is de- 
scribed by the density field 



n+{R) ^ yl(i?)exp 



A{R) exp— ^exp-^, 
(105) 

where is the first step the rotated k- vector is k+ = Rtfc, 
and in the second step we separated the fast oscillating 
factor exp{ik^ R/e^^'^) in the spirit of the multiscale ex- 
pansion. The first two factors are therefore the amplitude 
with respect to the basis set of the original k-vectors. 
The first important property of the box operator ia^Si^ 



□ exp 



0. 



(106) 



It implies that a pure crystal is a solution of the ampli- 
tude equations for arbitrary orientation. We note that 
we use both the notation of a scalar product (denoted by 
a dot •) and a matrix product (no multiplication symbol), 
i.e. a ■ b = d^b. From the definition of the operator we 
obtain 



□ exp 



iit'fMR 
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fct(R-I)fc+ ifct(R_I)(-Rt 



I)fc 



rl/2 



where we used the normalization condition = 1 and 
the orthogonality of the rotation matrix, R^R = 1. 

Next, we consider situations in which the whole solid is 
rotated, but the lattice orientation is kept in its original 
state. This is visualized in Fig. and d. We introduce 
rotated lattice vectors k- = R/c; notice that in compari- 
son to fc_|_ this vector is rotated in the opposite direction. 
Correspondingly, we define a rotated box operator 



1/2 



(107) 



If A{R) is a density wave amplitude, then ACRR) de- 
scribes the crystal with rotated shape, but the same lat- 
tice orientation. We obtain then 

□A(Ri?) = D-Ajj^^, (108) 

which expresses the equivalence of active and passive ro- 
tations: The rotated grain with original lattice orienta- 
tion (Fig. HJi) has the same properties as the original 



□A(Ri?) = A:1'R^ 



A 

RR 2 \RR 



RR ■ 



From the definition of the box operator follows imme- 
diately the product rule 



□ (/5) = fDg + gOf - ie^'\\I f) ■ (Vg) 



(109) 



We can also define a box operator with opposite rota- 
tion 



□+ = A:+ • V - 



le 
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(110) 



and from the product rule Eq. (|109p and Eq. ([T 
obtain the operator rotation rule 

- / - i^MRX ik'^MR ■ , 

□ I ^exp — ZY/2~ I ^ '^-'^P — TiT^u^A, 



we 
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(111) 



or, with the inverse rotation 



-/-f ik^M.'^R\ i^yVR- J 

°Me^P gi/2 I =exp— ^D.A. 



Analogous to Eq. (llOSp we have 

U+A{B.R)^UA\^^. 



(112) 



(113) 



Using Eqs. (|llip and (|113p we finally get the coordinate 
transformation rules towards a rotated frame of reference, 
i.e. the transition from a to e in Fig. 21 



□ |^^(Ri?) exp ^ 
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(114) 



and 



|^A(Ri?) exp 



ik^M.R 

ei/2 



exp ■ 



ik^yiR 

rl/2 



U'A\ 



RR ■ 



(115) 

The latter equation expresses the rotational invariance 
of the amplitude equations. Namely, if A{R) is a valid 
amplitude field, the field A(R,R) describes a rotated ma- 
terial, but still with the same lattice orientation as be- 
fore. The multiplication with the exponential factor cor- 
responds to the rotation of the lattice only, and therefore 
the expression in brackets on the left hand side describes 
the rotated material. As expressed by this relation, we 
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1.006 



1.004 



1.002 



Rotated interface 
Rotated lattice 
Rotated interface, no box operator 




TABLE 11: Solid-melt interfacial free energy for different in- 
terface orientations of bee iron. The value are given in i/rr? . 
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FIG. 5: Solid-liquid interface energy as function of orientation 
for the two-dimensional hexagonal system. For £213 = 0.1 the 
influence of the higher order term in the box operator gives 
only a small correction. The curves show results from a rota- 
tion of the interface normal vector, whereas the points show 
results for differently rotated lattices. The results coincide 
if the box operator is included; without it, a rotated lattice 
structure has a higher bulk energy and is therefore not in 
equilibrium with the liquid at the nominal melting tempera- 
ture. 7°; is the solid-liquid interfacial energy for 9 = Q and 
£213 = 0, i.e. without the higher order correction of the box 
operator. 



obtain the physical equivalence of the state, and therefore 
the rotational invariance. Notice that all local terms also 
acquire the same exponential rotation factor, and there- 
fore this factor cancels in the end in the homogeneous 
amplitude equations. 

We therefore conclude that also the free energy remains 
unchanged by a rotation of both the lattice and the mi- 
crostructure. This, in contrast, is not true without the 
corrective term of the box operator, since the operator 
fc • V is not rotational invariant'^ -''^^ . If we consider e.g. 
a planar solid-melt interface, which is stable exactly at 
the melting point, independent of the interface normal 
direction. The rotational invariant formulation with the 
box operator preserves this if the whole system is rotated 
(transition a to d in Fig.U]), where the amplitudes acquire 
apart from the rotation of the microstructure profile also 
the beats, i.e. A{R) A{IIR) exp(ifc'l'M^/ei/2). ^^j^g gj^_ 
ergy of the solid phase remains unchanged (this is trivial 
for the melt, since there the amplitudes vanish and are 
therefore always invariant). Without the higher order 
term in the box operator, however, the energy density 
of the solid increases spuriously, and therefore the solid 
would start to melt. 

To make this more transparent, we calculated 
the anisotropic surface energy density for the two- 
dimensional hexagonal system as function of the interface 
orientation, as shown in Fig. [5j The curves are obtained 
from equilibration runs to minimize the free energy with 
fixed k- vectors, but with different interface normal vec- 



Orientation 


without box operator 


with box operator 


100 


0.14414 


0.14392 


110 


0.14067 


0.14051 


111 


0.13576 


0.13643 



tors, as done in Refs. |42||46| . Since we assume a straight 
interface, all amplitudes depend only on the normal di- 
rection, and the problem becomes one-dimensional. The 
dashed curve shows the result without correction term in 
the box operator (i.e. formally setting e = 0), the solid 
curve for finite, small e. Both curves differ only very lit- 
tle, in agreement with the fact that the higher order term 
in the box operator gives only a small correction. Notice 
that since we do not rotate the reciprocal lattice vectors 
but only the normal vector, the solid-liquid interface re- 
mains stable at T = Tm- The points, in contrast, show 
data with rotated reciprocal lattice vectors and fixed in- 
terface normal. As expected, the results fall exactly onto 
the curve with fixed reciprocal lattice vectors and rotated 
interface normal. Without the box operator corrections, 
the equilibration would lead to a pure liquid (no phase 
coexistence), since the solid bulk energy would be raised 
artificially, thus making the solid unfavorable at the nom- 
inal melting temperature; this behavior would obviously 
be unphysical. We note that the simulations with ro- 
tated lattice vectors require the solution of the full two- 
dimensional problem, since the amplitudes depend now 
on both coordinates due to the beats of the exponential 
factor. 

Table |lT] lists the equilibrium interfacial free energies 
between solid and melt for bcc iron, using the parameters 
shown in tablelH Here we clearly see that the higher order 
term gives only a small correction to the values calculated 
in Ref. 42. 

It follows from the frame invariance of the crystal-melt 
interfacial free-energies, that the amplitude equation ap- 
proach should describe well the solidification of a poly- 
crystalline material from an undercooled melt. We illus- 
trate this here for the case of two-dimensional hexagonal 
crystals (see Fig. HI) . Several spherical seed crystals with 
different orientation are implanted into the melt phase 
and grow (provided that they exceed the critical radius). 
When the crystals meet, they form grain boundaries, 
which can consist of isolated dislocations for low angle 
grain boundaries or show a rather diffuse interface re- 
gion, which can be partially premelted. Notice that the 
defect distribution is not static but slowly evolves, since 
the dislocations interact with each other via long-range 
elastic forces. 



19 




FIG. 6: Illustration of polycrystalline isothermal solidification of two-dimensional, pure hexagonal crystals, modeled by the 
amplitude equations of the phase field model, as given in AppendixlX] The actual system is much larger than the magnification 
shown here. When interfaces between misoriented grains form, they generate isolated dislocations for subgrain boundaries and 
diffuse and premelted interfaces for high misorientations. The parameters are T = —0.002 and e2c ~ 0.1. 



C. Grain boundary energies and premelting 

Up to this point, the amplitude equations reflect the 
rotation invariance of the physical system correctly, and 
this is related to the fact that the melt is fully rotational 
invariant, since all amplitudes vanish there. The situ- 
ation becomes more complex if we consider a polycrys- 
tal. Let us consider a grain boundary, where the energy 
depends on the orientation of both crystals. E.g. for a 
hexagonal crystal, it is obvious that apart from the con- 
tinuous symmetries to which we paid attention so far. 



also discrete symmetries are important: If we rotate one 
of the adjacent crystals by 60°, it is in the same state 
again, and therefore the grain boundary energy has not 
changed - it exhibits a sixfold symmetry. 

The dependence of the grain boundary energy on mis- 
orientation is shown in Fig. [7] for a symmetric tilt bound- 
ary in a hexagonal crystal. The temperature is apprecia- 
bly below the melting point, in order to "stabilize" the 
grain boundary and to prevent a large separation of the 
grains due to premelting, which will be briefly discussed 
below. Starting from a dense-packed configuration (see 
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FIG. 7: Grain boundary energy as a function of misorienta- 
tion for symmetrical tilt boundaries in a hexagonal crystal for 
two different inclinations, normalized to the solid-liquid inter- 
facial energy. The dimensionless undercooling is T — —0.01, 
and €20 = 0.1. The results do not obey sixfold symmetry. 




e [degree] 

FIG. 8: Grain boundary energy as function of the misorienta- 
tion for symmetric tilt grain boundaries in bcc iron. The two 
curves correspond to [100] and [110] interface normals. The 
temperature is Tm — T — 80K. 
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FIG. 9: Grain boundary energy as a function of misorienta- 
tion for symmetrical tilt boundaries in a hexagonal crystal 
exactly at the melting point {T = Tm), normalized to twice 
the solid-liquid interfacial energy. The amplitude equation re- 
sults are in very good quantitative agreement with the PFC 
results from Ref . |23 for the same parameter £20 ~ 0.1. The 
boundary premelts for 6 larger than 9c ~ 10°. 



0° and 30° for the (j) = 0° inclination and 30° and 60° 
for the other cf) = 30° inclination (and similarly for bcc 
on either side of 45°.) A more detailed analysis for the 
simplest case of a smectic crystal is given in Appendix IB] 
Fig. [5] shows the grain boundary energy for small mis- 
orientations at a symmetrical grain boundary for the 
dense-packed crystal surfaces. In contrast to Fig. [71 
the temperature here is equal to the melting tempera- 
ture. Above a critical misorientation 6c ~ 10°, where 
7gh ~ 27s/, the grains premelt, and the thickness of the 
melt layer diverges logarithmically as the melting point 
is approached from below. The obtained data coincides 
well with the PFC simulations and a Read-Shockley fit, 
where the dislocation core radius is the only adjustable 
parameter. A more detailed investigation of grain bound- 
ary premelting in the context of the ampltiude equations 
will be discussed elsewhere. 



left panel of Fig. [Tj inclination cj) = 0) the misorienta- 
tion is increased, and we see that the grain boundary en- 
ergy increases monotonically. It therefore does not reflect 
the proper sixfold symmetry which would imply that 7gf, 
goes to zero for 6 = 60°. Conversely, starting from the 
(j) = 30° incliniation, the grain boundary does not "heal" 
if the dense-packed configuration is reached. A similar 
behavior is observed for bcc iron, where the amplitude 
equations do not obey the correct cubic symmetry, see 
Fig. ISl For the reasons explained in section I, the am- 
plitude equations are strictly valid only in the limit of 
small misorientations. However, for both hexagonal and 
bcc crystals the predictions remain approximately valid 
over roughly half the complete range allowed by the full 
crystal symmetry, e.g. is approximately valid between 



D. Shear-induced grain boundary coupling and 
sliding 

If a bicrystal is sheared in the direction parallel to a 
grain boundary, it migrates in a direction normal to the 
grain boundary plane for low temperatures^*^^^ , and this 
effect is contained in the amplitude equation formulation. 
In Fig. [ini the right crystal is sheared downwards. Mo- 
tion of the grain boundary by one lattice unit takes place 
during the time that an atom of the sheared right crystal 
needs to move until it matches the lattice of the left grain. 
Then the grain boundary shifts in normal direction with 
velocity v± = '(;s/[2 tan(6'/2)] for —tt/6 < 9 < tt/6, where 
Vg is the sliding velocity. This confirms the geometrical 
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model of coupling in Refs. 1681 - 1701 here applied to the 
hexagonal crystal symmetry. The results in Fig. [TT] also 
show that the amplitude equation approach reproduces 
the transition from coupled motion to sliding'''', where 
the latter is favored close enough to the melting point. 
A more detailed study of this transition in an amplitude 
equation framework is currently in progress. 

For the simulation we use a real space implementation, 
since we do not have periodic boundary conditions in 
the direction of the grain boundary normal. At the left 
boundary a; = we keep all amplitudes fixed in time 
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with M; = R; — I and the two-dimensional rotation ma- 
trix R; for the rotation of the left grain, compare also 
to equation ()103p and the notation introduced in Section 
IVIBI The boundary conditions for the right grain at 
X = X involve not only a rotation with R^ = R| in the 
opposite direction, but also a time-dependent displace- 
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where the displacement vector u has components = 
and Uy = 2e^ytX, with the strain rate e°j^ (defined 
here on the "slow" scale). Notice in particular that the 
elastic deformation factor involves the rotated principal 
reciprocal vectors R^^'^-'^ 



Acknowledgments 

This work was supported by DOE through grant DE- 
FG02-07ER46400 and the Computational Materials Sci- 
ence Network program. R.S. also acknowledges finan- 
cial support for the later part of this work of the Ger- 
man DFG grant SPP 1296 and from the industrial spon- 
sors of ICAMS, ThyssenKrupp Steel AG, Salzgitter Man- 
nesmann Forschung GmbH, Robert Bosch GmbH, Bayer 
Materials Science AG, Bayer Technology Services GmbH, 
Benteler AG and the state of North-Rhine- Westphalia. 



• • • . • 

• • • . ■ 

• a - • • 





• 



■000 

• • . • • . • 




FIG. 10: Shear-induced coupled motion of two grains. The 
left grain is fixed and the right one is slowly pulled downwards. 
The blue line shows the location of the grain boundary, the red 
and green lines the crystallographic planes. The atom which 
is marked by the green circle reaches the atomic plane of the 
left crystal between (b) and (c) and attaches to the left crys- 
tal (red circle). Through this mechanism, the grain boundary 
moves perpendicular to the pulling direction. Parameters are 
(in the notation of Appendix IX)) £213 =0.1, Ke2D = 1, with 
equal kinetic coefficients in the solid and the melt. The real 
space discretization is Ax = 0.25, timestep /S.t = 0.1, and the 
misorientation between the two crystals at the symmetric tilt 
is 23.28°. The shear rate is e^y = — 10~^, and the dimension- 
less undercooling T = —0.1. The snapeshots are taken from 
a to d at times 0.1, 9.0, 12.5 and 51.4, respectively. 
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FIG. 11: Sliding of two grains. The left grain is fixed, the 
right one is slowly pulled downwards. The red and green 
mark lines of the crystallographic planes. First, the system 
quickly equilibrates to a finite separation of the grains. Notice 
that the shear stress that builds up while the two grains are 
still connected also favors melting. Afterwards, the separa- 
tion and the interface positions do not change any more and 
only the right grain slides downward. The parameters are the 
same as in Fig. 1101 with the exception of the much smaller 
dimensionless undercooling T = —0.00012. 



Appendix A: Phase field crystal 

For completeness, we derive in this appendix the am- 
plitude equations from the phase field crystal model for 
two-dimensional hexagonal crystals using the same NWS 
type multiscale expansion as for bco^. 

We start from the dimensionless free energy functional 



dr( |[-62D + (V2 + l)2]^-,-iv/ ) (Al) 



The parameter €20 plays the same role as the scale sepa- 
ration parameter e as defined in section IIIIl Equilibrium 
requires that the chemical potential 
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'e2DV^+(V2-|-l)V + V'^ 
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is spatially constant. The density "0 is constant in the 
hquid, tp — tpi, thus we get the free energy density 
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We use a one-mode approximation for the solid, 
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'0(r) = f/js -I- Aj exp{ik^^^ ■ r) 
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with the following set of iV = 6 normalized principal 
reciprocal lattice vectors: 




The k- vectors used in the expansion above are a multiple 
of these vectors, and their length is determined by a free 
energy minimization below; in fact, we will obtain then 
fcO) = kU) . The amplitudes are real with 
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where the factor q is introduced to find the correct length 
of the reciprocal lattice vectors, which corresponds to the 
atomic spacing. We can then calculate the free energy 
(per unit cell) and minimize it with respect to q and A^. 
This gives^^ 
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where the ± sign is for positive and negative ips respec- 
tively; in accordance with Ref. |2^ we pick the nega- 
tive branch. Furthermore, we obtain q = \/3/2 (thus 
= fc(j)). Then the average free energy density in the 
sohd is 
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Coexistence between sohd and hquid demands the equal- 
ity of the chemical potentials, = ^g/i = dfs/i/dil)s/i^ 
and the grand potentials, w = ujs/i — Js/i ~ "^sjiP^^ of the 
two phases. 

In the spirit of a multiscale expansion we write the 
average densities as 



■f/'sO^D + V'sie2D + '0s2e2D + ■ ■ • (^10) 

i'iv>^2D + V';ie2D + i>i2^-iD + — (All) 



Up to the order t2D the chemical potential difference is 



which implies 
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From the next order term of the chemical potential bal- 
ance, ©(ej^^), we get a relation between %l^i2 and ■(/'s2, 
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-125V^,2v'-36i^2 

-90y^-36V,^o + 15^30 + 1200V',^o " 100 
+138.^-36V^2^ + 15V4o - 23G4V4o) 
/ (-125^-36V^2^ + 15). 
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Up to the order e^^^ the difference between the grand 
potentials vanishes, and from the order e\jj we get the 
result 
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(there is also another solution, l/-\/3, which we drop, 
since we concentrate here on the negative branch). Sim- 
ilarly, from the chemical potential balance at order e\jj 
we obtain tpn = tpsi = 0. 

Beyond the thermodynamical analysis above, which 
deals only with the spatially averaged quantities, we con- 
sider now explicitly the spatial oscillations of the density. 
The key is the separation of slow variables, denoted by 
capital letters R, and fast variables, r. They are related 



-* 1/2 -* 

by the expansion parameter £20, R = ^2d^- This trans- 
lates also to gradients, where we introduce two gradients, 
Vf and V^, where the first operator acts only on fast, 

the second on slow variables, thus V 
We therefore obtain the transformation rule 
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(V2-+l)2 + 4e^/,2(v2_+l)V,-.V^ 



-2e 



2D 



(V|-hl)V| + 2(V,.-V^) 
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where we skipped terms of order e^j^ and higher. 
We expand the field according to 

^(f) = MMd + M^^2D + Mr)4D +■■■■ (A16) 



Using the expansion for the averaged densities Eqs. (lAlOp 
and (jAlip . the multiscale phase equilibrium version of 
the PFC equation is 



-e2DV'+(V2 + l)V+?/'^ = AQ4^D+{'^l2~^sO+i^sQ)4D+- ' ' ' 

which is of course the same as for hcc^, but with a dif- 
ferent value Ipso- 
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At order we obtain the equation 
which is solved by 

N 

^o(r) = V.o + E^°(^)e''"''''- 

At order e2Di we have 

(V|+ 1)2^1 (r) = 0, 

with the solution 
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At order £2^^ we have for the first time also gradients 
with respect to the slow variables, 

(yj+l)^^j2 - V^o + 4(V,- V^) + ^0 = ^12 - ^sO + ^sO- 

(A21) 

All secular terms proportional to exp(ifc'-^' • r) must bal- 
ance each other, thus we get e.g. for the prefactor in 
front of exp(ifc(^) • r) 

- 4(fc(i) . W^fAl + 6V'.o^r^r + SV-fo^? 
+6[A°yl°* + v4°A°*]A° + 3A°A°*A^l = 0. 

For a pure solid, we can find solution with constant real 
amplitudes A°, with Af = A^ ^ yl° -4/V^, and 
A2 = —A'^; this solution minimizes the free energy given 
below with F = 0. 
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The above equation (and its counterpart for the pref- 
actors of the other exponential terms) can be obtained 
variationally from the following free energy: 

r N/2 



N/2 



+6'4;so{A'i*A°2*A°*+AlAlAl) 

/jV/2 \^ N/2 >, 

+3(EA,%-j (A22) 

where we treat real and imaginary of the complex fields 
A° as independent for the variation. The prefactor of 
the functional, -F^^D' ^'^^ determined by the equilib- 
rium condition, and will be determined by consideration 
of elastic deformations below. The summation up to N/2 
means that we sum over 1, 2, 3. Notice that the same ex- 
pression holds for bcc up to quadratic order, apart from 
the fact the = 12 there for the different set of prin- 
cipal reciprocal lattice vectors. The cubic and quartic 
terms all satisfy the conditions that only vectors, which 
form closed polygons, contribute; this corresponds to the 
appearance of a 5-function when the fast oscillations are 
integrated out. 

The generalization to a rotational invariant form is via 
the replacement 



9 R 



(A23) 



The equilibrium conditions, that lead to the cancella- 
tion of the secular terms above is 
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We can formulate a dynamical form of these equations 
by 



dt 
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with a kinetic coefficient K. Notice that the free energy 
decays monotonically with these evolution equations for 
the complex fields, and we therefore finally reach an equi- 
librium state. 

Analogous to the general expression for the elastic con- 
stants ([5^ we obtain 
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3 for two distinct pairs of indices 
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We note that these expressions are defined on the "fast 
scale" r. As expected, this case corresponds to isotropic 
elasticity, and the usual Lame coefficient and shear mod- 
ulus are 
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This corresponds to a (three-dimensional) Poisson ratio 
of z/ = 1/4. 

Finally, this allows to determine the energy scale 
by calculation of the elastic energy of a deformed solid, 
where all amplitudes are equal to In the energy 

density of the phase field crystal model, e.g. a (small) 
strain txx leads to an increase of the free energy den- 
sity by fs — fi — (24/185)e^j.e2Li to lowest order in £20, 
if all other strain components vanish. On the other 
hand, the same free energy density change (|A22p for 
the amplitude equations, written of the fast scale, is 
fei — '^ij^ij/'^ = (A + fi/2)e1x with the elastic con- 
stants given in Eqs. (IA27|) and (IA28p . Thus we have 
fei = (24/185)^20^62^, and therefore 



^2D 



<^2D- 



(A29) 



A deviation from the melting temperature and the cou- 
pling to thermodynamic alloy models is achieved via an 
additional free energy term 



dR (t>T 
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with a dimensionless temperature deviation T from the 
melting point; the interpolating "phase field" is defined 
as in Eqs. ^ and (|35|) . 

We can eliminate the phase field crystal parameters by 
rescaling the equations using 



Aj — J^]/As 
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and introduction of another (small) dimensionless param- 
eter 



1 



^2D 



3 i^so ~ 1 



^2D- 



(A32) 



Similarly, the length scales are scaled with this new pa- 



rameter, X = t^j^x and 
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Then the free energy becomes 
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Appendix B: Breakdown of rotational invariance 

To illustrate the breakdown of rotational invariance, 
we consider the simple case of a smectic crystal, which 
is described by only two antiparallel principal reciprocal 
lattice vectors and only one complex amplitude u. A pure 
crystal is then described by the amplitude u{X) = 1 and 
the density variation 



5n{x) — ueyi'p{ix) + u* exp(— ix) 



(Bl) 



with X = If the crystal is rotated by 180°, 

i.e. u{X) = exp(— 2i2;), it recovers its original state ac- 
cording to Eq. (jB2|) with the same density 5n. Due to 
the high rotation angle, the spacing between the beats in 
the amplitude is half the lattice spacing, thus a fine dis- 
cretization is necessary for a numerical implementation. 
Notice that both states are purely one-dimensional. 

This system is described through a (dimensionless) 
one-dimensional free energy functional 
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, □ 

ul")" and /t — L{T — TM)h{\u\'^)/TM ■ In equilibrium, 
a one-dimensional solution is therefore described by the 
ordinary differential equation 
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where the prime ' denotes differentiation with respect to 
X. Obviously, both above one-dimensional states satisfy 
this equation, therefore the rotational invariance of the 
bulk states holds. 

Although both amplitudes describe the same density, a 
defect free interface cannot be formed between them. In 
particular, for T = Tm a "grain boundary" between the 
state u = exp(— 2ix) for x —oo and u = 1 for a: — >■ oo 
would premelt, and a melt layer forms between the two 
"grains" . In an undercooled state, T < Tm, Fig.[T2]shows 
the reconstructed equilibrium density, which seems to be 
defect-free. In the same figure, also the corresponding 
amplitude and the free energy density (without the con- 
tribution from the thermal tilt) are shown. At the inter- 
face, the phase &{X) of the amplitude, u = \u\ exp(i9), 
changes smoothly, see Fig. [131 Obviously, a finite inter- 
face forms between the two grains, and it is accompanied 
by a finite energy density, thus a spurious finite grain 
boundary energy is found. 




FIG. 12: Density variation 5n(X) and real part of the am- 
plitude u[X) of a smectic crystal that consists of one "grain" 
that is not rotated [X > 0) and one "grain" that is rotated 
by 180° [X < 0). The density seems to indicate a perfectly 
healed crystal, but nevertheless an interface energy is asso- 
ciated with the interface, leading to a finite energy density 
/*: + fp in the interface region and therefore a nonvanishing 
grain boundary energy (solid black area). The parameters 
used here are e = 0.1 and L{T - Tm)/Tm = -0.1. 
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FIG. 13: Phase and modulus of the amplitude u in the inter- 
face region, corresponding to the results and parameters in 
Fig.fll 



Appendix C: Periodic boundary conditions 

For purposes of numerical modeling using spectral 
methods^, periodic boundary conditions are advanta- 
geous. This becomes restrictive here, since we need the 
periodicity for each order parameter. For the simplest 
case of a liquid or a solid that is not rotated with re- 
spect to the "natural" orientation of the set of principal 
reciprocal lattice vectors, the amplitudes are constant in 
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each phase, and therefore no periodicity constraints arise. 
This is also true for coexistence of sohd and hquid. 

The situation becomes more complex if rotated crystals 
are involved, because then the density wave amplitudes 
acquire a periodic modulation. We look at the 2D hexag- 
onal system with size X xY first and use the notation 
of Appendix |^ Let us consider the situation of a single 
solid phase, then the amplitudes of the rotated crystal 
are given by 



(CI) 



Periodicity requires therefore at the horizontal boundary 
at a; = and x = X 



exp 




exp 




(C2) 



for j = I . . . N and arbitrary vertical coordinate y. 
Hence, 
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= 27rn,- 



(C3) 



with integer numbers Uj. Summation over the first 
N/2 = 3 principal reciprocal lattice vectors (see 
Eq. (|A5p ). which form a closed triangle, therefore gives 



ni + n2 + 713 = 0. 



(C4) 



(Since the fields that are associated to the other prin- 
cipal reciprocal lattice vectors are complex conjugate to 
the previous, their periodicity does not give additional 
conditions). From the definition of the rotation matrix 
Eq. (|104p we therefore get the two conditions 
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From these two equations we obtain immediately 
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V3(cos6»- 1) = — (2ni +712). (C7) 
X 



Subsequent division by Eq. (jC6[) and use of trigonometric 
identities yields 



r- 6 2711 + 712 
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2 772 ■ 
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which defines a discrete set of admissible of lattice rota- 
tion angles 0, given through the integer number ni,772. 
The corresponding system length is then 

27rn2 



X = 



sin 9 ' 
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which obviously becomes large for low angle rotations. 
The analogous expressions for periodicity in y direction 
in a system of height Y are 
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where 7711,7772 are also integer numbers, the system 
height is then 



Y = 



27r77i2 

COS 6* - 1 ' 



(Cll) 



Now, both conditions (jCSP and (jClOP must be satisfied, 
therefore giving 
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for which integer solutions have to be found. Since we 
are interested in finding small system sizes, i.e. small val- 
ues of 772 and 7772, we can choose 7772 = — 1 and 7711 < 0. 
Using for example 771 — —mi — 27772 and 772 = 2777i + 7772 
therefore gives a solution of the above equation with the 
angle given by Eq. (ICIOI) . and the system sizes follow 
from Eqs. (IC9p and (jClip . For specific cases, also solu- 
tions with smaller systems sizes can be found. 

Although the requirement of periodicity for a bcc crys- 
tal seems to be more stringent at a first glance due to the 
higher number of order parameters, this complexity is 
significantly reduced by the ability to form different sets 
of closed polygons of principal reciprocal lattice vectors. 
First, we have now six conditions of the type (ICSp . but 
from the formation of closed triangles we get the integer 
relations 



'^011 + "101 - "110 = 0) 
-77011 + "110 - "loi = 0> 
-77oii + 77110 - "101 = 0, 
0, 



'oil 



noi 



'110 



(CIS) 
(C14) 
(C15) 
(C16) 



but only three of them are independent. Also, the condi- 
tions of closed quadrilaterals does not provide additional 
independent information. 

The periodicity condition for [Oil] is explicitly in anal- 
ogy to Eq. jCS]) 
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sin OX — 27r77oii 



(C17) 



and similar for [101] and [110 
1 



V2 



{cos9-l)X = 27r77ioi, (C18) 



^(cose*- 1)X - ^sin6lX = 27r77iio, (C19) 
V 2 V 2 

from which get the additional (independent) integer re- 
lation 



"no - "101 - "-on = 0. 



(C20) 



Therefore, only two integer numbers can be chosen inde- 
pendently, e.g. 77101 and 77011, and we finally arrive at the 
expressions 
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Similarly, periodicity in y direction implies the conditions 
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with the only two independent integer numbers toioi and 
moil- The equality of the angles according to Eqs. (|C21[) 
and (|C23[) demands therefore integer solutions of the 
equation 
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If we assume translation invariance of the amplitudes 
in z direction (i.e. no periodicity condition) , the problem 
becomes effectively two-dimensional. Notice, however, 
that the reconstructed density waves still have a periodic 
modulation in that direction, and also the atoms are not 
bound to stay in the xy plane. In fact, all displacements 
Ux{x,y), Uy{x,y) and Uz{x,y), that depend only on in- 
plane coordinates, can be described by a two-dimensional 
formulation of the amplitude equations. 



For a discussion concerning periodicity in systems with 
a grain boundary we refer to Ref. 12^ 
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